Construct Single-Hierarchical P/NBD Model for Long Timeframe Synthetic Data

In this workbook we construct the non-hierarchical P/NBD models on the synthetic data with the longer timeframe.

1 Load and Construct Datasets

1.1 Load Long-Timeframe Synthetic Transaction Data

We now want to load the CD-NOW transaction data.

Code

customer_cohortdata_tbl <- read_rds("data/longsynth_customer_cohort_data_tbl.rds")
customer_cohortdata_tbl |> glimpse()

Rows: 5,000
Columns: 5
$ customer_id     <chr> "LFC201001_0001", "LFC201001_0002", "LFC201001_0003", …
$ cohort_qtr      <chr> "2010 Q1", "2010 Q1", "2010 Q1", "2010 Q1", "2010 Q1",…
$ cohort_ym       <chr> "2010 01", "2010 01", "2010 01", "2010 01", "2010 01",…
$ first_tnx_date  <dttm> 2010-01-01 16:16:30, 2010-01-03 16:32:54, 2010-01-04 …
$ total_tnx_count <int> 12, 1, 1, 2, 2, 3, 2, 7, 5, 5, 4, 8, 1, 1, 1, 31, 1, 1…

Code

customer_transactions_tbl <- read_rds("data/longsynth_transaction_data_tbl.rds")
customer_transactions_tbl |> glimpse()

Rows: 42,244
Columns: 7
$ customer_id   <fct> LFC201001_0001, LFC201001_0002, LFC201001_0003, LFC20100…
$ tnx_timestamp <dttm> 2010-01-01 16:16:30, 2010-01-03 16:32:54, 2010-01-04 22…
$ tnx_dow       <fct> Fri, Sun, Mon, Tue, Wed, Thu, Sun, Sun, Mon, Mon, Tue, W…
$ tnx_month     <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, J…
$ tnx_week      <chr> "00", "00", "01", "01", "01", "01", "01", "01", "02", "0…
$ invoice_id    <chr> "T20100101-0001", "T20100103-0001", "T20100104-0001", "T…
$ tnx_amount    <dbl> 160.87, 1.30, 96.11, 813.07, 85.58, 35.07, 12.23, 70.48,…

Code

customer_subset_id <- read_rds("data/longsynth_customer_subset_ids.rds")
customer_subset_id |> glimpse()

 Factor w/ 5000 levels "LFC201001_0001",..: 2 8 10 14 16 17 27 30 33 35 ...

We re-produce the visualisation of the transaction times we used in previous workbooks.

Code

plot_tbl <- customer_transactions_tbl |>
  group_nest(customer_id, .key = "cust_data") |>
  filter(map_int(cust_data, nrow) > 3) |>
  slice_sample(n = 30) |>
  unnest(cust_data)

ggplot(plot_tbl, aes(x = tnx_timestamp, y = customer_id)) +
  geom_line() +
  geom_point() +
  labs(
      x = "Date",
      y = "Customer ID",
      title = "Visualisation of Customer Transaction Times"
    ) +
  theme(axis.text.y = element_text(size = 10))

1.2 Load Derived Data

Code

obs_fitdata_tbl   <- read_rds("data/longsynth_obs_fitdata_tbl.rds")
obs_validdata_tbl <- read_rds("data/longsynth_obs_validdata_tbl.rds")

customer_fit_stats_tbl <- obs_fitdata_tbl |>
  rename(x = tnx_count)

1.3 Load Subset Data

We also want to construct our data subsets for the purposes of speeding up our valuations.

Code

customer_fit_subset_tbl <- obs_fitdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_fit_subset_tbl |> glimpse()

Rows: 1,000
Columns: 6
$ customer_id    <fct> LFC201001_0002, LFC201001_0008, LFC201001_0010, LFC2010…
$ first_tnx_date <dttm> 2010-01-03 16:32:54, 2010-01-10 08:26:52, 2010-01-11 2…
$ last_tnx_date  <dttm> 2010-01-03 16:32:54, 2010-06-24 02:01:46, 2010-02-28 0…
$ tnx_count      <dbl> 0, 6, 4, 30, 0, 0, 4, 2, 2, 11, 0, 6, 3, 1, 1, 1, 0, 4,…
$ t_x            <dbl> 0.0000000, 23.5332238, 6.7394202, 30.6992466, 0.0000000…
$ T_cal          <dbl> 625.7586, 624.8069, 624.5765, 624.2573, 624.1894, 624.0…

Code

customer_valid_subset_tbl <- obs_validdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_valid_subset_tbl |> glimpse()

Rows: 1,000
Columns: 3
$ customer_id       <fct> LFC201001_0002, LFC201001_0008, LFC201001_0010, LFC2…
$ tnx_count         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ tnx_last_interval <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

We now use these datasets to set the start and end dates for our various validation methods.

Code

dates_lst <- read_rds("data/longsynth_simulation_dates.rds")

use_fit_start_date <- dates_lst$use_fit_start_date
use_fit_end_date   <- dates_lst$use_fit_end_date

use_valid_start_date <- dates_lst$use_valid_start_date
use_valid_end_date   <- dates_lst$use_valid_end_date

We now split out the transaction data into fit and validation datasets.

Code

customer_fit_transactions_tbl <- customer_transactions_tbl |>
  filter(
    customer_id %in% customer_subset_id,
    tnx_timestamp >= use_fit_start_date,
    tnx_timestamp <= use_fit_end_date
    )
  
customer_fit_transactions_tbl |> glimpse()

Rows: 9,521
Columns: 7
$ customer_id   <fct> LFC201001_0002, LFC201001_0008, LFC201001_0010, LFC20100…
$ tnx_timestamp <dttm> 2010-01-03 16:32:54, 2010-01-10 08:26:52, 2010-01-11 23…
$ tnx_dow       <fct> Sun, Sun, Mon, Thu, Thu, Fri, Thu, Sat, Mon, Tue, Tue, F…
$ tnx_month     <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, J…
$ tnx_week      <chr> "00", "01", "02", "02", "02", "02", "03", "03", "04", "0…
$ invoice_id    <chr> "T20100103-0001", "T20100110-0002", "T20100111-0002", "T…
$ tnx_amount    <dbl> 1.30, 70.48, 830.32, 29.04, 33.90, 45.65, 36.21, 78.96, …

Code

customer_valid_transactions_tbl <- customer_transactions_tbl |>
  filter(
    customer_id %in% customer_subset_id,
    tnx_timestamp >= use_valid_start_date,
    tnx_timestamp <= use_valid_end_date
    )
  
customer_valid_transactions_tbl |> glimpse()

Rows: 738
Columns: 7
$ customer_id   <fct> LFC202105_0029, LFC201908_0005, LFC201906_0021, LFC20211…
$ tnx_timestamp <dttm> 2022-01-01 06:56:30, 2022-01-01 17:47:30, 2022-01-02 01…
$ tnx_dow       <fct> Sat, Sat, Sun, Sun, Sun, Sun, Mon, Mon, Mon, Mon, Mon, M…
$ tnx_month     <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, J…
$ tnx_week      <chr> "00", "00", "00", "00", "00", "00", "01", "01", "01", "0…
$ invoice_id    <chr> "T20220101-0004", "T20220101-0008", "T20220102-0001", "T…
$ tnx_amount    <dbl> 328.72, 213.08, 14.54, 17.96, 67.87, 19.41, 2.05, 11.22,…

Finally, we want to extract the first transaction for each customer, so we can add this data to assess our models.

Code

customer_initial_tnx_tbl <- customer_fit_transactions_tbl |>
  slice_min(n = 1, order_by = tnx_timestamp, by = customer_id)

customer_initial_tnx_tbl |> glimpse()

Rows: 1,000
Columns: 7
$ customer_id   <fct> LFC201001_0002, LFC201001_0008, LFC201001_0010, LFC20100…
$ tnx_timestamp <dttm> 2010-01-03 16:32:54, 2010-01-10 08:26:52, 2010-01-11 23…
$ tnx_dow       <fct> Sun, Sun, Mon, Thu, Thu, Fri, Thu, Sat, Tue, Sat, Sat, W…
$ tnx_month     <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, F…
$ tnx_week      <chr> "00", "01", "02", "02", "02", "02", "03", "03", "04", "0…
$ invoice_id    <chr> "T20100103-0001", "T20100110-0002", "T20100111-0002", "T…
$ tnx_amount    <dbl> 1.30, 70.48, 830.32, 29.04, 33.90, 45.65, 36.21, 78.96, …

We now expand out these initial transactions so that we can append them to our simulations.

Code

sim_init_tbl <- customer_initial_tnx_tbl |>
  transmute(
    customer_id,
    draw_id       = list(1:n_sim),
    tnx_timestamp,
    tnx_amount
    ) |>
  unnest(draw_id)

sim_init_tbl |> glimpse()

Rows: 2,000,000
Columns: 4
$ customer_id   <fct> LFC201001_0002, LFC201001_0002, LFC201001_0002, LFC20100…
$ draw_id       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
$ tnx_timestamp <dttm> 2010-01-03 16:32:54, 2010-01-03 16:32:54, 2010-01-03 16…
$ tnx_amount    <dbl> 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1…

Before we start on that, we set a few parameters for the workbook to organise our Stan code.

Code

stan_modeldir <- "stan_models"
stan_codedir  <-   "stan_code"

2 Fit First Hierarchical Lambda Model

Our first hierarchical model puts a hierarchical prior around the mean of our population \(\lambda\) - lambda_mn.

Once again we use a Gamma prior for it.

2.1 Compile and Fit Stan Model

We now compile this model using CmdStanR.

Code

pnbd_onehierlambda_stanmodel <- cmdstan_model(
  "stan_code/pnbd_onehier_lambda.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )

We then use this compiled model with our data to produce a fit of the data.

Code

stan_modelname <- "pnbd_long_onehierlambda1"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.25,
    hier_lambda_mn_p2 = 1,

    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 1.00,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_long_onehierlambda1_stanfit <- pnbd_onehierlambda_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_long_onehierlambda1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_long_onehierlambda1_stanfit <- read_rds(stanfit_object_file)
}

pnbd_long_onehierlambda1_stanfit$print()

  variable      mean    median    sd   mad        q5       q95 rhat ess_bulk
 lp__      -97358.53 -97359.00 74.33 76.50 -97482.11 -97239.46 1.00      716
 lambda_mn      0.26      0.26  0.01  0.01      0.25      0.26 1.00     1304
 lambda[1]      0.04      0.03  0.01  0.01      0.02      0.05 1.00     2550
 lambda[2]      0.15      0.09  0.19  0.10      0.00      0.50 1.00     1930
 lambda[3]      0.15      0.09  0.17  0.10      0.00      0.49 1.01     1265
 lambda[4]      0.24      0.18  0.19  0.15      0.03      0.64 1.00     1876
 lambda[5]      0.26      0.21  0.20  0.17      0.04      0.66 1.00     2190
 lambda[6]      0.19      0.17  0.12  0.11      0.05      0.42 1.00     2405
 lambda[7]      0.20      0.16  0.16  0.12      0.03      0.52 1.00     2562
 lambda[8]      0.23      0.21  0.09  0.08      0.11      0.40 1.00     2087
 ess_tail
     1311
     1569
     1235
      979
      881
      973
     1350
     1164
     1174
     1578

 # showing 10 of 13728 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

We have some basic HMC-based validity statistics we can check.

Code

pnbd_long_onehierlambda1_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda1-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda1-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda1-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda1-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

2.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "lambda_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_long_onehierlambda1_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_long_onehierlambda1_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

Finally, we want to check out the energy diagnostic, which is often indicative of problems with the posterior mixing.

Code

pnbd_long_onehierlambda1_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)

2.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_long_onehierlambda1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_long_onehierlambda1_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_long_onehierlambda1",
  fit_end_dttm        = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm    = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm      = use_valid_end_date   |> as.POSIXct(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 3110
  )

pnbd_long_onehierlambda1_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_long_onehierlambda1_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_long_onehierlambda1_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_long_onehierlambda1_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_long_onehierlambda1_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_long_onehierlambda1_assess_valid_simstats_tbl.rds"

2.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_long_onehierlambda1_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

2.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_long_onehierlambda1_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

As for our short time frame data, overall our model is working well.

3 Fit Second Hierarchical Lambda Model

In this model, we are going with a broadly similar model but we are instead using a different mean for our hierarchical prior.

3.1 Fit Stan Model

We now want to fit the model to our data using this alternative prior for lambda_mn.

Code

stan_modelname <- "pnbd_long_onehierlambda2"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.50,
    hier_lambda_mn_p2 = 1,

    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 1.00,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_long_onehierlambda2_stanfit <- pnbd_onehierlambda_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_long_onehierlambda2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_long_onehierlambda2_stanfit <- read_rds(stanfit_object_file)
}

pnbd_long_onehierlambda2_stanfit$print()

  variable      mean    median    sd   mad        q5       q95 rhat ess_bulk
 lp__      -97363.25 -97363.50 71.55 69.46 -97481.34 -97244.88 1.00      526
 lambda_mn      0.26      0.26  0.00  0.01      0.25      0.26 1.00     1360
 lambda[1]      0.03      0.03  0.01  0.01      0.02      0.05 1.00     2295
 lambda[2]      0.14      0.09  0.17  0.09      0.01      0.51 1.00     1726
 lambda[3]      0.15      0.09  0.18  0.10      0.00      0.52 1.00     1725
 lambda[4]      0.24      0.20  0.19  0.16      0.03      0.63 1.00     2067
 lambda[5]      0.26      0.21  0.21  0.17      0.04      0.66 1.00     1981
 lambda[6]      0.19      0.17  0.12  0.11      0.05      0.43 1.00     1961
 lambda[7]      0.20      0.16  0.16  0.13      0.03      0.51 1.00     1746
 lambda[8]      0.23      0.22  0.09  0.09      0.10      0.40 1.00     2580
 ess_tail
     1101
     1661
     1443
      830
     1011
     1111
     1262
     1422
     1051
     1495

 # showing 10 of 13728 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

We have some basic HMC-based validity statistics we can check.

Code

pnbd_long_onehierlambda2_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda2-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda2-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda2-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehierlambda2-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

3.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "lambda_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_long_onehierlambda2_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_long_onehierlambda2_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

Finally, we want to check out the energy diagnostic, which is often indicative of problems with the posterior mixing.

Code

pnbd_long_onehierlambda2_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)

3.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_long_onehierlambda2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_long_onehierlambda2_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_long_onehierlambda2",
  fit_end_dttm        = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm    = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm      = use_valid_end_date   |> as.POSIXct(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 3120
  )

pnbd_long_onehierlambda2_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_long_onehierlambda2_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_long_onehierlambda2_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_long_onehierlambda2_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_long_onehierlambda2_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_long_onehierlambda2_assess_valid_simstats_tbl.rds"

3.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_long_onehierlambda2_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

3.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_long_onehierlambda2_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

As for our short time frame data, overall our model is working well.

4 Fit First Hierarchical Mu Model

We now construct the same hierarchical model but based around mu_mn.

4.1 Compile and Fit Stan Model

We compile this model using CmdStanR.

Code

pnbd_onehiermu_stanmodel <- cmdstan_model(
  "stan_code/pnbd_onehier_mu.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )

We then use this compiled model with our data to produce a fit of the data.

Code

stan_modelname <- "pnbd_long_onehiermu1"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_mu_mn_p1 = 0.50,
    hier_mu_mn_p2 = 1.00,

    lambda_mn     = 0.25,
    lambda_cv     = 1.00,

    mu_cv         = 1.00

    )

if(!file_exists(stanfit_object_file)) {
  pnbd_long_onehiermu1_stanfit <- pnbd_onehiermu_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_long_onehiermu1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_long_onehiermu1_stanfit <- read_rds(stanfit_object_file)
}

pnbd_long_onehiermu1_stanfit$print()

  variable      mean    median    sd   mad        q5       q95 rhat ess_bulk
 lp__      -93164.34 -93163.75 73.79 74.65 -93284.83 -93042.68 1.00      732
 mu_mn          0.11      0.11  0.00  0.00      0.10      0.11 1.00      527
 lambda[1]      0.04      0.03  0.01  0.01      0.02      0.05 1.01     2415
 lambda[2]      0.15      0.08  0.19  0.10      0.01      0.51 1.00     1998
 lambda[3]      0.15      0.08  0.19  0.10      0.01      0.52 1.00     1707
 lambda[4]      0.24      0.19  0.19  0.15      0.03      0.58 1.00     2715
 lambda[5]      0.25      0.20  0.20  0.16      0.04      0.64 1.00     1889
 lambda[6]      0.19      0.17  0.12  0.10      0.05      0.43 1.00     2048
 lambda[7]      0.20      0.16  0.17  0.12      0.03      0.51 1.00     1516
 lambda[8]      0.23      0.21  0.09  0.08      0.11      0.39 1.00     2405
 ess_tail
     1160
     1221
     1312
     1056
      918
     1322
     1363
     1126
     1044
     1508

 # showing 10 of 13728 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

We have some basic HMC-based validity statistics we can check.

Code

pnbd_long_onehiermu1_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu1-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu1-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu1-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu1-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

4.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "mu_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_long_onehiermu1_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_long_onehiermu1_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

Finally, we want to check out the energy diagnostic, which is often indicative of problems with the posterior mixing.

Code

pnbd_long_onehiermu1_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)

4.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_long_onehiermu1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_long_onehiermu1_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_long_onehiermu1",
  fit_end_dttm        = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm    = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm      = use_valid_end_date   |> as.POSIXct(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 3130
  )

pnbd_long_onehiermu1_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_long_onehiermu1_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_long_onehiermu1_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_long_onehiermu1_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_long_onehiermu1_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_long_onehiermu1_assess_valid_simstats_tbl.rds"

4.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_long_onehiermu1_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

4.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_long_onehiermu1_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

As for our short time frame data, overall our model is working well.

5 Fit Second Hierarchical Lambda Model

In this model, we are going with a broadly similar model but we are instead using a different mean for our hierarchical prior.

5.1 Fit Stan Model

We now want to fit the model to our data using this alternative prior for lambda_mn.

Code

stan_modelname <- "pnbd_long_onehiermu2"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_mu_mn_p1 = 0.25,
    hier_mu_mn_p2 = 1.00,

    lambda_mn     = 0.25,
    lambda_cv     = 1.00,

    mu_cv         = 1.00
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_long_onehiermu2_stanfit <- pnbd_onehiermu_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_long_onehiermu2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_long_onehiermu2_stanfit <- read_rds(stanfit_object_file)
}

pnbd_long_onehiermu2_stanfit$print()

  variable      mean    median    sd   mad        q5       q95 rhat ess_bulk
 lp__      -93168.62 -93167.00 74.49 76.06 -93293.11 -93047.89 1.01      577
 mu_mn          0.11      0.11  0.00  0.00      0.10      0.11 1.01      532
 lambda[1]      0.03      0.03  0.01  0.01      0.02      0.05 1.00     2406
 lambda[2]      0.14      0.09  0.16  0.09      0.01      0.44 1.00     2292
 lambda[3]      0.14      0.08  0.16  0.09      0.00      0.44 1.00     1722
 lambda[4]      0.24      0.18  0.20  0.15      0.03      0.63 1.00     2028
 lambda[5]      0.25      0.20  0.20  0.16      0.04      0.64 1.00     2325
 lambda[6]      0.19      0.17  0.11  0.10      0.05      0.40 1.00     1967
 lambda[7]      0.20      0.16  0.16  0.13      0.03      0.51 1.00     1786
 lambda[8]      0.22      0.21  0.09  0.09      0.10      0.38 1.00     2608
 ess_tail
     1005
     1005
     1375
     1516
     1077
     1356
     1323
     1125
     1003
     1496

 # showing 10 of 13728 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

We have some basic HMC-based validity statistics we can check.

Code

pnbd_long_onehiermu2_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu2-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu2-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu2-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_long_onehiermu2-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

5.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "mu_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_long_onehiermu2_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_long_onehiermu2_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

Finally, we want to check out the energy diagnostic, which is often indicative of problems with the posterior mixing.

Code

pnbd_long_onehiermu2_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)

5.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_long_onehiermu2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_long_onehiermu2_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_long_onehiermu2",
  fit_end_dttm        = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm    = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm      = use_valid_end_date   |> as.POSIXct(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 3140
  )

pnbd_long_onehiermu2_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_long_onehiermu2_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_long_onehiermu2_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_long_onehiermu2_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_long_onehiermu2_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_long_onehiermu2_assess_valid_simstats_tbl.rds"

5.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_long_onehiermu2_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

5.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_long_onehiermu2_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)

$total_plot

$quant_plot

As for our short time frame data, overall our model is working well.

6 Compare Model Outputs

We have looked at each of the models individually, but it is also worth looking at each of the models as a group.

Code

calculate_simulation_statistics <- function(file_rds) {
  simdata_tbl <- read_rds(file_rds)
  
  multicount_cust_tbl <- simdata_tbl |>
    filter(sim_tnx_count > 0) |>
    count(draw_id, name = "multicust_count")
  
  totaltnx_data_tbl <- simdata_tbl |>
    count(draw_id, wt = sim_tnx_count, name = "simtnx_count")
  
  simstats_tbl <- multicount_cust_tbl |>
    inner_join(totaltnx_data_tbl, by = "draw_id")
  
  return(simstats_tbl)
}

Code

obs_fit_customer_count <- customer_fit_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_valid_customer_count <- customer_valid_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_fit_total_count <- customer_fit_subset_tbl |>
  pull(tnx_count) |>
  sum()

obs_valid_total_count <- customer_valid_subset_tbl |>
  pull(tnx_count) |>
  sum()


obs_stats_tbl <- tribble(
  ~assess_type, ~name,               ~obs_value,
  "fit",        "multicust_count",   obs_fit_customer_count,
  "fit",        "simtnx_count",      obs_fit_total_count,
  "valid",      "multicust_count",   obs_valid_customer_count,
  "valid",      "simtnx_count",      obs_valid_total_count
  )


model_assess_tbl <- dir_ls("data", regexp = "pnbd_long_(one|fixed).*_assess_.*simstats") |>
  enframe(name = NULL, value = "file_path") |>
  filter(str_detect(file_path, "_assess_model_", negate = TRUE)) |>
  mutate(
    model_label = str_replace(file_path, "data/pnbd_long_(.*?)_assess_.*", "\\1"),
    assess_type = if_else(str_detect(file_path, "_assess_fit_"), "fit", "valid"),
    
    sim_data = map(
      file_path, calculate_simulation_statistics,
      
      .progress = "calculate_simulation_statistics"
      )
    )

model_assess_tbl |> glimpse()

Rows: 14
Columns: 4
$ file_path   <fs::path> "data/pnbd_long_fixed1_assess_fit_simstats_tbl.rds", …
$ model_label <chr> "fixed1", "fixed1", "fixed2", "fixed2", "fixed3", "fixed3"…
$ assess_type <chr> "fit", "valid", "fit", "valid", "fit", "valid", "fit", "va…
$ sim_data    <list> [<tbl_df[2000 x 3]>], [<tbl_df[2000 x 3]>], [<tbl_df[2000…

Code

model_assess_summstat_tbl <- model_assess_tbl |>
  select(model_label, assess_type, sim_data) |>
  unnest(sim_data) |>
  pivot_longer(
    cols = !c(model_label, assess_type, draw_id)
    ) |>
  group_by(model_label, assess_type, name) |>
  summarise(
    .groups = "drop",
    
    mean_val = mean(value),
    p10 = quantile(value, 0.10),
    p25 = quantile(value, 0.25),
    p50 = quantile(value, 0.50),
    p75 = quantile(value, 0.75),
    p90 = quantile(value, 0.90)
    )

model_assess_summstat_tbl |> glimpse()

Rows: 28
Columns: 9
$ model_label <chr> "fixed1", "fixed1", "fixed1", "fixed1", "fixed2", "fixed2"…
$ assess_type <chr> "fit", "fit", "valid", "valid", "fit", "fit", "valid", "va…
$ name        <chr> "multicust_count", "simtnx_count", "multicust_count", "sim…
$ mean_val    <dbl> 639.2355, 8441.8860, 58.1700, 733.4730, 701.0205, 4670.491…
$ p10         <dbl> 621.0, 7672.9, 53.0, 641.0, 683.9, 4299.9, 41.0, 367.0, 60…
$ p25         <dbl> 630.00, 7995.75, 56.00, 684.75, 692.00, 4462.00, 43.00, 40…
$ p50         <dbl> 639.0, 8425.5, 58.0, 733.5, 701.0, 4664.0, 45.0, 445.0, 62…
$ p75         <dbl> 649.00, 8859.00, 61.00, 783.00, 710.00, 4867.25, 48.00, 48…
$ p90         <dbl> 656.0, 9243.0, 63.0, 826.0, 720.0, 5071.0, 50.0, 518.0, 63…

Code

#! echo: TRUE

ggplot(model_assess_summstat_tbl) +
  geom_errorbar(
    aes(x = model_label, ymin = p10, ymax = p90), width = 0
    ) +
  geom_errorbar(
    aes(x = model_label, ymin = p25, ymax = p75), width = 0, linewidth = 3
    ) +
  geom_hline(
    aes(yintercept = obs_value),
    data = obs_stats_tbl, colour = "red"
    ) +
  scale_y_continuous(labels = label_comma()) +
  expand_limits(y = 0) +
  facet_wrap(
    vars(assess_type, name), scale = "free_y"
    ) +
  labs(
    x = "Model",
    y = "Count",
    title = "Comparison Plot for the Different Models"
    ) +
  theme(
    axis.text.x = element_text(angle = 20, vjust = 0.5, size = 8)
    )

6.1 Write Assessment Data to Disk

We now want to save the assessment data to disk.

Code

model_assess_tbl |> write_rds("data/assess_data_pnbd_long_onehier_tbl.rds")

R Environment

Code

options(width = 120L)
sessioninfo::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16)
 os       Ubuntu 22.04.3 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Dublin
 date     2023-11-20
 pandoc   3.1.1 @ /usr/local/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version    date (UTC) lib source
 abind            1.4-5      2016-07-21 [1] RSPM (R 4.3.0)
 arrayhelpers     1.1-0      2020-02-04 [1] RSPM (R 4.3.0)
 backports        1.4.1      2021-12-13 [1] RSPM (R 4.3.0)
 base64enc        0.1-3      2015-07-28 [1] RSPM (R 4.3.0)
 bayesplot      * 1.10.0     2022-11-16 [1] RSPM (R 4.3.0)
 bit              4.0.5      2022-11-15 [1] RSPM (R 4.3.0)
 bit64            4.0.5      2020-08-30 [1] RSPM (R 4.3.0)
 bridgesampling   1.1-2      2021-04-16 [1] RSPM (R 4.3.0)
 brms           * 2.20.4     2023-09-25 [1] RSPM (R 4.3.0)
 Brobdingnag      1.2-9      2022-10-19 [1] RSPM (R 4.3.0)
 cachem           1.0.8      2023-05-01 [1] RSPM (R 4.3.0)
 callr            3.7.3      2022-11-02 [1] RSPM (R 4.3.0)
 checkmate        2.3.0      2023-10-25 [1] RSPM (R 4.3.0)
 cli              3.6.1      2023-03-23 [1] RSPM (R 4.3.0)
 cmdstanr       * 0.6.0.9000 2023-11-07 [1] Github (stan-dev/cmdstanr@a13c798)
 coda             0.19-4     2020-09-30 [1] RSPM (R 4.3.0)
 codetools        0.2-19     2023-02-01 [2] CRAN (R 4.3.1)
 colorspace       2.1-0      2023-01-23 [1] RSPM (R 4.3.0)
 colourpicker     1.3.0      2023-08-21 [1] RSPM (R 4.3.0)
 conflicted     * 1.2.0      2023-02-01 [1] RSPM (R 4.3.0)
 cowplot        * 1.1.1      2020-12-30 [1] RSPM (R 4.3.0)
 crayon           1.5.2      2022-09-29 [1] RSPM (R 4.3.0)
 crosstalk        1.2.0      2021-11-04 [1] RSPM (R 4.3.0)
 curl             5.1.0      2023-10-02 [1] RSPM (R 4.3.0)
 digest           0.6.33     2023-07-07 [1] RSPM (R 4.3.0)
 directlabels   * 2023.8.25  2023-09-01 [1] RSPM (R 4.3.0)
 distributional   0.3.2      2023-03-22 [1] RSPM (R 4.3.0)
 dplyr          * 1.1.3      2023-09-03 [1] RSPM (R 4.3.0)
 DT               0.30       2023-10-05 [1] RSPM (R 4.3.0)
 dygraphs         1.1.1.6    2018-07-11 [1] RSPM (R 4.3.0)
 ellipsis         0.3.2      2021-04-29 [1] RSPM (R 4.3.0)
 evaluate         0.22       2023-09-29 [1] RSPM (R 4.3.0)
 fansi            1.0.5      2023-10-08 [1] RSPM (R 4.3.0)
 farver           2.1.1      2022-07-06 [1] RSPM (R 4.3.0)
 fastmap          1.1.1      2023-02-24 [1] RSPM (R 4.3.0)
 forcats        * 1.0.0      2023-01-29 [1] RSPM (R 4.3.0)
 fs             * 1.6.3      2023-07-20 [1] RSPM (R 4.3.0)
 furrr          * 0.3.1      2022-08-15 [1] RSPM (R 4.3.0)
 future         * 1.33.0     2023-07-01 [1] RSPM (R 4.3.0)
 generics         0.1.3      2022-07-05 [1] RSPM (R 4.3.0)
 ggdist           3.3.0      2023-05-13 [1] RSPM (R 4.3.0)
 ggplot2        * 3.4.4      2023-10-12 [1] RSPM (R 4.3.0)
 globals          0.16.2     2022-11-21 [1] RSPM (R 4.3.0)
 glue           * 1.6.2      2022-02-24 [1] RSPM (R 4.3.0)
 gridExtra        2.3        2017-09-09 [1] RSPM (R 4.3.0)
 gtable           0.3.4      2023-08-21 [1] RSPM (R 4.3.0)
 gtools           3.9.4      2022-11-27 [1] RSPM (R 4.3.0)
 hms              1.1.3      2023-03-21 [1] RSPM (R 4.3.0)
 htmltools        0.5.6.1    2023-10-06 [1] RSPM (R 4.3.0)
 htmlwidgets      1.6.2      2023-03-17 [1] RSPM (R 4.3.0)
 httpuv           1.6.12     2023-10-23 [1] RSPM (R 4.3.0)
 igraph           1.5.1      2023-08-10 [1] RSPM (R 4.3.0)
 inline           0.3.19     2021-05-31 [1] RSPM (R 4.3.0)
 jsonlite         1.8.7      2023-06-29 [1] RSPM (R 4.3.0)
 knitr            1.44       2023-09-11 [1] RSPM (R 4.3.0)
 labeling         0.4.3      2023-08-29 [1] RSPM (R 4.3.0)
 later            1.3.1      2023-05-02 [1] RSPM (R 4.3.0)
 lattice          0.21-8     2023-04-05 [2] CRAN (R 4.3.1)
 lifecycle        1.0.3      2022-10-07 [1] RSPM (R 4.3.0)
 listenv          0.9.0      2022-12-16 [1] RSPM (R 4.3.0)
 loo              2.6.0      2023-03-31 [1] RSPM (R 4.3.0)
 lubridate      * 1.9.3      2023-09-27 [1] RSPM (R 4.3.0)
 magrittr       * 2.0.3      2022-03-30 [1] RSPM (R 4.3.0)
 markdown         1.11       2023-10-19 [1] RSPM (R 4.3.0)
 Matrix           1.5-4.1    2023-05-18 [2] CRAN (R 4.3.1)
 matrixStats      1.0.0      2023-06-02 [1] RSPM (R 4.3.0)
 memoise          2.0.1      2021-11-26 [1] RSPM (R 4.3.0)
 mime             0.12       2021-09-28 [1] RSPM (R 4.3.0)
 miniUI           0.1.1.1    2018-05-18 [1] RSPM (R 4.3.0)
 munsell          0.5.0      2018-06-12 [1] RSPM (R 4.3.0)
 mvtnorm          1.2-3      2023-08-25 [1] RSPM (R 4.3.0)
 nlme             3.1-162    2023-01-31 [2] CRAN (R 4.3.1)
 parallelly       1.36.0     2023-05-26 [1] RSPM (R 4.3.0)
 pillar           1.9.0      2023-03-22 [1] RSPM (R 4.3.0)
 pkgbuild         1.4.2      2023-06-26 [1] RSPM (R 4.3.0)
 pkgconfig        2.0.3      2019-09-22 [1] RSPM (R 4.3.0)
 plyr             1.8.9      2023-10-02 [1] RSPM (R 4.3.0)
 posterior      * 1.4.1      2023-03-14 [1] RSPM (R 4.3.0)
 prettyunits      1.2.0      2023-09-24 [1] RSPM (R 4.3.0)
 processx         3.8.2      2023-06-30 [1] RSPM (R 4.3.0)
 promises         1.2.1      2023-08-10 [1] RSPM (R 4.3.0)
 ps               1.7.5      2023-04-18 [1] RSPM (R 4.3.0)
 purrr          * 1.0.2      2023-08-10 [1] RSPM (R 4.3.0)
 quadprog         1.5-8      2019-11-20 [1] RSPM (R 4.3.0)
 QuickJSR         1.0.7      2023-10-15 [1] RSPM (R 4.3.0)
 R6               2.5.1      2021-08-19 [1] RSPM (R 4.3.0)
 Rcpp           * 1.0.11     2023-07-06 [1] RSPM (R 4.3.0)
 RcppParallel     5.1.7      2023-02-27 [1] RSPM (R 4.3.0)
 readr          * 2.1.4      2023-02-10 [1] RSPM (R 4.3.0)
 reshape2         1.4.4      2020-04-09 [1] RSPM (R 4.3.0)
 rlang          * 1.1.1      2023-04-28 [1] RSPM (R 4.3.0)
 rmarkdown        2.25       2023-09-18 [1] RSPM (R 4.3.0)
 rstan            2.32.3     2023-10-15 [1] RSPM (R 4.3.0)
 rstantools       2.3.1.1    2023-07-18 [1] RSPM (R 4.3.0)
 rstudioapi       0.15.0     2023-07-07 [1] RSPM (R 4.3.0)
 rsyslog        * 1.0.3      2023-05-08 [1] RSPM (R 4.3.0)
 scales         * 1.2.1      2022-08-20 [1] RSPM (R 4.3.0)
 sessioninfo      1.2.2      2021-12-06 [1] RSPM (R 4.3.0)
 shiny            1.7.5.1    2023-10-14 [1] RSPM (R 4.3.0)
 shinyjs          2.1.0      2021-12-23 [1] RSPM (R 4.3.0)
 shinystan        2.6.0      2022-03-03 [1] RSPM (R 4.3.0)
 shinythemes      1.2.0      2021-01-25 [1] RSPM (R 4.3.0)
 StanHeaders      2.26.28    2023-09-07 [1] RSPM (R 4.3.0)
 stringi          1.7.12     2023-01-11 [1] RSPM (R 4.3.0)
 stringr        * 1.5.0      2022-12-02 [1] RSPM (R 4.3.0)
 svUnit           1.0.6      2021-04-19 [1] RSPM (R 4.3.0)
 tensorA          0.36.2     2020-11-19 [1] RSPM (R 4.3.0)
 threejs          0.3.3      2020-01-21 [1] RSPM (R 4.3.0)
 tibble         * 3.2.1      2023-03-20 [1] RSPM (R 4.3.0)
 tidybayes      * 3.0.6      2023-08-12 [1] RSPM (R 4.3.0)
 tidyr          * 1.3.0      2023-01-24 [1] RSPM (R 4.3.0)
 tidyselect       1.2.0      2022-10-10 [1] RSPM (R 4.3.0)
 tidyverse      * 2.0.0      2023-02-22 [1] RSPM (R 4.3.0)
 timechange       0.2.0      2023-01-11 [1] RSPM (R 4.3.0)
 tzdb             0.4.0      2023-05-12 [1] RSPM (R 4.3.0)
 utf8             1.2.4      2023-10-22 [1] RSPM (R 4.3.0)
 V8               4.4.0      2023-10-09 [1] RSPM (R 4.3.0)
 vctrs            0.6.4      2023-10-12 [1] RSPM (R 4.3.0)
 vroom            1.6.4      2023-10-02 [1] RSPM (R 4.3.0)
 withr            2.5.1      2023-09-26 [1] RSPM (R 4.3.0)
 xfun             0.40       2023-08-09 [1] RSPM (R 4.3.0)
 xtable           1.8-4      2019-04-21 [1] RSPM (R 4.3.0)
 xts              0.13.1     2023-04-16 [1] RSPM (R 4.3.0)
 yaml             2.3.7      2023-01-23 [1] RSPM (R 4.3.0)
 zoo              1.8-12     2023-04-13 [1] RSPM (R 4.3.0)

 [1] /usr/local/lib/R/site-library
 [2] /usr/local/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Code

options(width = 80L)